//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Psi_Cell_Field.i.hh
 * \author Jeremy Roberts
 * \date   Aug 17, 2011
 * \brief  Psi_Cell_Field inline member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//


#ifndef PSI_CELL_FIELD_I_HH_
#define PSI_CELL_FIELD_I_HH_

#include <cmath>

namespace slabtran
{

// Constructor
template<class Equations, class D>
Psi_Cell_Field<Equations,D>::Psi_Cell_Field(SP_mesh     mesh,
                                            int         number_angles,
                                            double      initial_value)
  :   Cell_Field<double> (mesh,
                          number_angles * Equations::number_unknowns(),
                          initial_value)
    , d_number_angles(number_angles)
{
  // checks?
}

// ----------------------------------------------------------------------------
// Mutable accessors

// 3d psi accessor via operator(); 2d and 1d calls this.
template<class Equations, class D>
inline typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int j, int k, int unknown, int angle)
{
  Require(angle >= 0);
  Require(angle < d_number_angles);
  Require(unknown >= 0);
  Require(unknown < Equations::number_unknowns());

  // In sweeps, we access all the cell/unknown values for
  // a given angle.  That is for each angle, we sweep over
  // cell and then unknown within a cell.  The most cache-
  // efficient way to store then seems to be
  // angle-->cell-->unknown.  Note, this overrides the
  // default indexing of Cell_Field.
  int index = angle * d_mesh->number_cells() * Equations::number_unknowns() +
              d_mesh->convert(i, j, k) * Equations::number_unknowns() +
              unknown;
  // This cardinal index is essentially:
  //   number of angles before this one * number values per angle +
  //   number of cells before this one * number values per cell +
  //   index of the unknown in this cell.

  // We got an index; now make sure it's valid.
  Ensure(valid_index(index));
  return d_values[index];
}

// 2d interface via operator()
template<class Equations, class D>
inline typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int j, int unknown, int angle)
{
  return (*this)(i, j, 0, unknown, angle);
}

// 1d interface via operator()
template<class Equations, class D>
inline typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int unknown, int angle)
{
  return (*this)(i, 0, 0, unknown, angle);
}

// ----------------------------------------------------------------------------
// Const accessors

// 3d psi accessor via operator(); 2d and 1d calls this.
template<class Equations, class D>
inline const typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int j, int k, int unknown, int angle) const
{
  Require(angle >= 0);
  Require(angle < d_number_angles);
  Require(unknown >= 0);
  Require(unknown < Equations::number_unknowns());

  // In sweeps, we access all the cell/unknown values for
  // a given angle.  That is for each angle, we sweep over
  // cell and then unknown within a cell.  The most cache-
  // efficient way to store then seems to be
  // angle-->cell-->unknown.  Note, this overrides the
  // default indexing of Cell_Field.
  int index = angle * d_mesh->number_cells() * Equations::number_unknowns() +
              d_mesh->convert(i, j, k) * Equations::number_unknowns() +
              unknown;
  // This cardinal index is essentially:
  //   number of angles before this one * number values per angle +
  //   number of cells before this one * number values per cell +
  //   index of the unknown in this cell.

  // We got an index; now make sure it's valid.
  Ensure(valid_index(index));
  return d_values[index];
}

// 2d interface via operator()
template<class Equations, class D>
inline const typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int j, int unknown, int angle) const
{
  return (*this)(i, j, 0, unknown, angle);
}

// 1d interface via operator()
template<class Equations, class D>
inline const typename Psi_Cell_Field<Equations,D>::value_type&
Psi_Cell_Field<Equations,D>::
operator()(int i, int unknown, int angle) const
{
  return (*this)(i, 0, 0, unknown, angle);
}

} // end namespace slabtran

#endif /* PSI_CELL_FIELD_I_HH_ */

//---------------------------------------------------------------------------//
//              end of Psi_Cell_Field.i.hh
//---------------------------------------------------------------------------//
